rm(list= ls())

library(tidyverse)
library(lfe)
library(broom)

## Define rank_fun

rank_fun <- function (x) {
  su = sort(unique(x))
  for (i in 1:length(su)) x[x == su[i]] = i
  return(x)
}

## Load Data and Define Outcomes
## Note: In BTW, right_total_party = AfD + CDU/CSU

bt <- readRDS('data/federal_elections_clean.RDS') %>%
  mutate(year = date.id) %>% 
  filter(year %in% c(2009, 2013, 2017)) %>% 
  group_by(ags) %>% 
  mutate(period = rank_fun(year)) %>% 
  filter(state %in% c('Nordrhein-Westfalen', 'Bayern'))

## Nr of pre and post treatment periods

id_treat <- bt %>% 
  filter(post == 1) %>% 
  filter(period == min(period)) %>% 
  pull(period) %>% 
  unique()

## Create Period*Treated Indicator Variables 

bt <- fastDummies::dummy_cols(bt, select_columns = 'period') %>% 
  mutate_at(vars(contains('period_')), ~ .*treated) 

## Generate Model Formulas 

outcomevars <- c('right_total_party',
                 'left_total_party')
periodvars <- names(bt)[str_detect(names(bt), 'period_')]

## Exclude last pre-treatment period 

periodvars <- periodvars[-which(periodvars == 'period_2')]

## Interact Period * State: allow for state-specific trends 

bt$state_period <- paste(bt$period, bt$state, sep = '')

## Run models

results <- lapply(outcomevars, function(dv){
  
  model_formula <- as.formula(paste0(dv, ' ~', paste0(periodvars, collapse = '+'), 
                                     '| state_period + ags  |0 | AA.district.id'))
  
  model_est <- felm(model_formula, data = bt) %>%
    broom::tidy(conf.int = T) %>% 
    filter(str_detect(term, 'period')) %>% 
    mutate(outcome = paste(dv))
  
}) %>% reduce(rbind) 

## Rename the period variables

results <- results %>% 
  mutate(period = parse_number(term) - (as.numeric(id_treat) - 1)) %>%
  filter(str_detect(outcome, 'total'))

## More renaming

results <- results %>% 
  mutate(outcome = recode(outcome, 
                          `left_total_party` = 'Total left vote',
                          `right_total_party_noafd` = 'Total right vote (excl. AfD)',
                          `right_total_party` = 'Total right vote',
                          `turnout_party` = 'Turnout')) 

## Add zero for the last pre-treatment period

zero_period <- data.frame(term = rep('period_2', 2),
                          estimate = rep(0, 2),
                          conf.low = rep(0, 2),
                          conf.high =rep(0, 2),
                          outcome = unique(results$outcome),
                          period = rep(0, 2))

## Append

results <- bind_rows(results, zero_period) %>%
  mutate(period = as.factor(period))

## Plot

myLoc <- 
  (which(levels(factor(results$period)) == "1") +
     which(levels(factor(results$period)) == "0")) / 
  2

#### Figure A.10: Leads and lags: federal elections ####

p1 <- ggplot(results, aes(x = factor(period), y = estimate, 
                          ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 0, linetype = 'dotted', color = 'grey40') +
  geom_vline(xintercept = myLoc, linetype = 'dotted', color = 'grey40') +
  geom_errorbar(width = 0) +
  geom_point(fill = 'white', shape = 21) +
  theme_bw() + 
  ylab('Coefficient Estimate\n(percentage points)') + 
  xlab('') + 
  facet_wrap(~ outcome) + 
  scale_color_grey(name = '') + 
  scale_x_discrete(breaks = c('-1', '0', '1'), 
                   labels = c('2009', '2013', '2017')) +
  ylim(c(-8, 6))
p1

## Clear and load full data set again 

rm(list = ls()[!ls() == 'rank_fun'])

bt <- readRDS('data/federal_elections_clean.RDS') %>%
  mutate(year = date.id) %>% 
  group_by(ags) %>% 
  mutate(period = rank_fun(year)) %>% 
  filter(state %in% c('Nordrhein-Westfalen', 'Bayern')) %>% 
  mutate(right_total_party_nofdp = cdu_csu_party + ifelse(!is.na(afd_party),
                                                          afd_party,
                                                          0))
outcomevars <- c('right_total_party',
                 'left_total_party')

#### Trend plot ####

## Note: In BTW, right_total_party = AfD + CDU/CSU

tp_df <- bt %>% 
  group_by(treated, year, state) %>% 
  summarise(across(all_of(outcomevars), \(x) mean(x, na.rm = TRUE))) %>% 
  ungroup() %>% 
  pivot_longer(cols = -c("year", "treated", "state"), values_to = "share",
               names_to = "outcome") %>% 
  mutate(outcome = recode(outcome, 
                          `right_total_party_nofdp` = 'Total right vote (excl. FDP)',
                          `cdu_csu_party` = 'Total right vote (excl. AfD and FDP)',
                          `left_total_party` = 'Total left vote',
                          `right_total_party_noafd` = 'Total right vote (excl. AfD)',
                          `right_total_party` = 'Total right vote',
                          `turnout_party` = 'Turnout')) %>% 
  filter(str_detect(outcome, 'Total')) %>% 
  mutate(state = ifelse(state == 'Bayern',
                        'Bavaria',
                        'North Rhine-Westphalia'))

## Plot

ggplot(tp_df %>% 
         filter(outcome %in% c("Total left vote",
                             "Total right vote")), 
       aes(year, share, factor(treated))) +
  geom_line(aes(linetype = factor(treated))) +
  geom_point(aes(shape = factor(treated)), fill = 'white') +
  geom_vline(xintercept = 2015, linetype = 'dotted') +
  facet_grid(outcome~state) +
  theme_bw() +
  scale_x_continuous(breaks = unique(bt$year)) +
  xlab("Election") +
  ylab("Vote share (%)") +
  scale_linetype_discrete(name = '', labels = c("Restricted access",
                                                "Liberalized access")) +
  scale_shape_manual(name = '', labels = c("Restricted access",
                                           "Liberalized access"),
                     values = c(21, 22)) +
  theme(legend.position = "bottom", legend.key.width = unit(1, 'cm'))
